Effect of bending rigidity on the knotting of a 
polymer under tension t 

Richard Matthews,*'* Ard A. Louis,^ and Christos N. Likos* 

Faculty of Physics, University of Vienna, Boltzmanngasse 5, A-1090 Vienna, Austria, and Rudolf 
Peierls Centre for Theoretical Physics, 1 Keble Road, Oxford 0X1 3NP, England 

E-mail: richard.matthews@univie.ac.at 
Abstract 

A coarse-grained computational model is used to investigate how the bending rigidity of a 
polymer under tension affects the formation of a trefoil knot. Thermodynamic integration tech- 
niques are applied to demonstrate that the free-energy cost of forming a knot has a minimum 
at non-zero bending rigidity. The position of the minimum exhibits a power-law dependence 
on the applied tension. For knotted polymers with non-uniform bending rigidity, the knots 
preferentially localize in the region with a bending rigidity that minimizes the free-energy. 

Type II topoisomerases are enzymes that may knot or unknot DNA by introducing a transient 
break in both strands of one DNA duplex and passing a second duplex through it. One of their 
key biological functions is to regulate the level of knotting in the genome. 1 Type II topoisomerases 
tend to act preferentially on certain sequences in DNA. 2 There is evidence that sites that are more 
frequently cleaved tend to be located in or next to parts of the genome called scaffold associated 
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regions or matrix attachment regions, ~ 5 which are typically several hundred base pairs long and 
rich in adenine (A) and thymine (T), two of the nucleotides in DNA. Further, a specific sequence 
evolved in vitro, which was preferentially cleaved by a certain type II topoisomerase, was highly 
AT-rich. 2 

It is believed that AT-rich sequences are more flexible than random ones. 5-8 For example, the 
work of Okonogi et al. 7 suggests that a sequence of AT repeats is about 20% more flexible than a 
random sequence. An earlier study suggested that such an AT rich sequence can have a persistence 
length less than half that of a GC rich sequence. 6 Scipioni et al. 8 used scanning force microscopy 
to observe a correlation between AT-rich parts of a DNA fragment and flexibility. Further, Masilah 
et al. 5 found that there is a preferentially large opening of the base-pairs immediately adjacent to 
a preferentially cleaved site. This opening was found to be dependent on the sequence context. 
Opening of base-pairs (bubble formation) can lead to greatly increased local flexibility. 9 Very high 
flexibility at the topoisomerase II cleavage sites is probably necessary because the enzyme enforces 
a large bend in DNA when it binds to it. 10 

An intriguing question arises as to whether the correlation between the positions of cleavage 
sites and DNA flexibility could be important in the regulation of knotting. For example, could 
the variation of bending stiffness help to localize knots near cleavage sites, thus expediting their 
removal? Here we make a first step towards understanding these issues by using a simple bead- 
spring polymer model to investigate how the free energy cost of forming a knot, kFknomng, varies 
with polymer bending stiffness and how this influences the position of a knot within a polymer 
of non-uniform flexibility. In this work, we simulate only the trefoil knot, 3i, but our general 
arguments do not depend on the particular topology. Previous work 12 on how the action of type 
II topoisomerase may be guided by bent geometries of DNA has been performed, but variable 
bending stiffness was not considered. 

The case of polymers under tension is biologically relevant because the action of enzymes 
during processes such as transcription applies forces to DNA. 12 ' 13 In general, for polymers in a 
good solvent with bending stiffness, A, under tension, x, there are three main contributions to 
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AFknotting'- the reduction in entropy due the self-confinement of the polymer in the knotted region; 
the increase in bending energy due to the curvature enforced by the knot; and the work done 
against the tension in reducing the extension of the polymer, necessary to give free length for knot 
formation. 

We consider how ^Fknotting varies with A for fixed x. We identify two length scales: that 
associated with the bending stiffness, Ia ~A/(£gr), and that associated with the size of the knotted 
region, 4„ 0/ (A), which depends on A. When Ia <C hnot(A) the main effect of increasing A will be 
to decrease the entropic cost of knotting and /^knotting will decrease with A. Previous work on 
fully flexible chains (A = 0), has found knots to be weakly localized, 14-16 N^not ~ N', where Nk no t 
is the number of monomers in the knot, N the total number in the polymer, and < t < 1 . 14 By 
applying scaling arguments based on the blob picture to interpret the results of simulations of 
polymers under tension, Farago et al., 14 estimated t = 0.4 ± 0. 1. A later study used two methods, 
including one based on closing subsections of the polymer and calculating a knot invariant, to find 
t ~ 0.75. 15 The discrepancy between the two estimates may be attributed to the relatively short 
polymers used in the earlier work. 15 Knot localization has been observed experimentally 17 but 
is found to disappear with confinement. 18 A free energy calculation for an open, linear polymer 
found no evidence of a metastable knot size. 19 

In the flexible regime, a polymer under tension will form a linear series of blobs of Nt, ~ 
(kgT /x) l / v monomers each, where V ~ 3/5. 20 The series of blobs cannot be knotted and so the 
knot resides within one blob. Treating this blob as an independent polymer, we expect l^not to be 
determined by the entropic localization of the knot and the number of monomers participating to 
the knot to scale, accordingly, as N^not ~ (kgT /x)^ v . By employing the simulation techniques 
and knot-identification algorithm to be presented shortly, we have determined the dependence of 
Nknot on x for a flexible polymer ofN = 256 beads of size o each. The results in Figure 1 indeed 
show a power-law dependence. By fitting to this data, we estimate that t = 0.43 ±0.01, which is 
consistent with the value found by Farago et al., 14 as expected given the relatively short chains 
used. Concomitantly, the knot size in fully flexible chains scales as l\ wot (0) ~ N^ not ~ (kgT fx) 1 . 
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Figure 1: Variation of the number of beads forming the knot, Nk no t with tension, x for N = 256 
bead flexible polymers. The solid line is a fit to the data with slope —0.71 ±0.01. Errorbars were 
estimated by performing three independent repeats of the simulations 

For Ia 3> I knot (A) the size of the knot will be dominated by the interplay of bending energy 
and tension and ^knotting will increase with A. We therefore expect a minimum of &Fknotting{A) 
at a value of A determined by Ia ~ I knot (A)- As the dependence of I knot (A) on T is not known, 
we replace I knot (A) with ^ oi (0) to find what the likely form of the dependence of the bending 
stiffness for which the free energy cost is minimal, A„„„, on x is. Using the results obtained above, 
a power-law dependence is obtained: 

Amin ~ X~'. (1) 

Of course, the replacement of hnot{A) with lk no t(ty m the relationship Ia ~ hnot(A) is an ap- 
proximation which is expected to break down precisely in the region of validity of this equality. 
On the other hand, a power-law dependence Nknot ~ ^ A is a reasonable assumption also for the 
case A^0, thus we anticipate a relationship of the form of Eq. (1) to hold also for A ^ 0, albeit 
with some exponent tA^t. 

For very large values of A, we expect the knot to form a single loop with the all crossings close 
to each other. 21 Assuming the thickness of the polymer is small compared to the loop, we expect 
AFknotting to be approximately given by 21 

^knotting = VSn 2 Ax. (2) 
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For lower A, the form of kFknottine ma y not be so easily deduced. At the crossover this is par- 
ticularly difficult because here we expect the bending length and self-confinement length to be 
approximately equal. For this case, a scaling form of the confinement free energy is not avail- 
able. 22 

We next study the consequences of these predictions with computer simulations. In what fol- 
lows, we first outline the technical details of our approach, we then present results on &Fknotting> 
before investigating the positional probability distribution of knots in polymers of non-uniform 
flexibility. We primarily simulate single chains ofN = 256 beads of size o in a simulation box of 
volume V = 2.048 x 10 5 a 3 with periodic boundaries: unless otherwise stated, all results are for 
these parameters. The polymers are connected to themselves across the periodic boundaries in the 
x-direction. A constant tension is simulated by including in the potential a term proportional to the 
x-length of the box, L x and allowing L x to vary. The advantage of this approach is that there are no 
free ends so that, as long as chain crossings are prevented, unknotting will never occur. 

The simulation of the polymer is carried through for the following interaction potential: 
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where r ;j = rj — r,, is the vector from bead i to bead j, located at position vectors r ( - and r 7 , 
respectively, whereas r,j denotes a unit vector. The first term sets the bending stiffness, which 
may be varied along the chain using the parameter Ki, giving a bending stiffness of A = KjO for 
the z'th bead. The second term applies a tension, T. The third and fourth terms are spring and 
excluded volume terms respectively, H is the Heaviside step function which truncates the Lennard- 
Jones potential to be purely repulsive. We choose £ = k^T ', k = 30kgT /o 2 and Rq = 1.5(7, which 
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prevents the chain from crossing itself and so conserves topology. 

We simulate using a Monte Carlo (MC) algorithm, 23 which comprises two types of moves. To 
simulate a given tension, moves that attempt to change L x , whilst rescaling L y and L z to keep V 
fixed and also applying a corresponding transformation to all particle coordinates, are included. 
Displacements of the polymer beads are made using the Hybrid MC method, 24 where trial states 
are generated using Molecular Dynamics (MD). During the MD trajectories, L x is fixed, the tension 
term is not included in the Hamiltonian used to calculate the forces. Collective motions of the 
polymer beads are more easily captured in this way than by local, single bead moves. 

To calculate kFknotting f° r a given tension, T, we simulate systems with all K\ set to the same 
value, K. We simulate two sets of systems, one with linear topology and one with knotted poly- 
mers. The systems within one set span a range of rigidities from k = up to the desired value. 
For each of those values, we calculate the average ^^0- By numerically integrating from 
K = 0, we obtain the relative free energy as a function of K", 23 AF a (?c) = F a (fc) — F a (0), where a 
stands for either 'knot' or 'linear'. To fully determine ^knotting we would need to perform an inte- 
gration between unknotted and knotted states. However, since we are interested in the relative cost 
of knotting for different bending stiffnesses, we simply calculate AFk not ti n g( K ) — AF / t„ offin g(0) = 
Aiw(K) -AFiinear(K) instead. 

To improve the efficiency of our calculation of AFk nott j ng (K) — AFk nott i ng (0) we implemented 
the most computationally intensive part of our simulation algorithm on a GPU using CUDA, which 
allows for a high degree of parallelism but is restrictive in terms of the homogeneity of the parallel 
calculations. 25 Whilst a standard local-move MC algorithm would be difficult to implement on a 
GPU, 25 the most time-consuming part of our algorithm is calculating the MD trajectories to pro- 
duce trial states for the Hybrid MC. The MD integration may be straightforwardly performed on a 
GPU. We simulate all systems for a given x and topology in parallel, performing force calculations 
and integration steps on the GPU. As a simple alternative to a cell-list we reduce the number of 
pair separations calculated by exploiting the connectivity of the polymer, which guarantees the 
maximum separation of two beads within a section: by comparing the center of mass positions of 
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two sections we can determine whether beads within them may interact. Random number genera- 
tion and other MC moves were performed on the CPU. To help reduce correlation times we added 
parallel tempering 23 swaps between systems with different K. 

(a)l 
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Figure 2: Schematic depiction of the knot-finding process, (a) The polymer is divided into sections 
by finding points along its contour - indicated by the dashed lines - at which there is a boundary 
between regions where only one strand crosses the y-z plane and those where multiple strands do. 
Regions in which there are multiple crossers are identified, these are indicated by the shaded areas. 
They may be closed and the Alexander polynomial calculated to identify which of them contains 
the knot, (b) - (d) Subsequently, a finer determination of the knot position may be achieved by 
taking the knot-containing section and considering subsections of it. These are closed by extending 
the polymer in the x-direction, as shown by the dotted lines. The Alexander polynomial may then 
be calculated for each of these. The section with the correct Alexander polynomial that contains the 
least number of beads is taken as containing the knot, (b) - (d) show a few examples of subsections. 
The subsection shown in (c) would be identified: that in (b) is contains more beads and that in (d) 
has the wrong polynomial. 



For simulations considering the positional probability distribution or size of the knot, it is 
necessary to determine the knotted section of the polymer. We applied a method, summarized 
in Figure 2, based on calculating the Alexander polynomial, 1 1 A^(x), at x = — 2 for polymer sub- 
sections. 26 Since the polymer is extended in the x-direction by the tension, there will usually be 
jc-positions at which only one part of polymer crosses the y-z-plane. Regions that are bounded by 
such points are considered. Only one will have the correct A^— 2). The more exact position is 
then found by taking subsections of this region, closing them with extensions in the zbc-direction, 
and finding the shortest with the correct A^(— 2). The center of this section is taken as the knot 
position and the number of beads it contains as the knot size. This is the same method applied for 
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the determination of N^ mt for flexible chains earlier in this paper. 

Our procedure may occasionally result in a false identification of a knot due to extra crossings 
included by the closing sections. However, in previous studies the rate of such errors was found 
to be low and to usually involve sections larger than truly knotted ones. 26 We thus do not expect 
such pitfalls to significantly affect our results but we refer the interested reader to an in-depth 
consideration of such schemes. 27 We also found that, occasionally, no jc-positions with only one 
crossing of the y-z-plane were found. In this case, the knot position was not identified and so 
these configurations were neglected. The rate of such configurations was < 1% for all the results 
presented. As a further check we verified that, for the knot size results, if instead of neglecting the 
configurations, a knot size equal to the total polymer size was added, the final averages were not 
changed by more than the errorbars. Simulations with knot-finding were performed with the same 
MC algorithm as for the free energy calculations. However, due to the computational cost of the 
knot-finding algorithm, which would be difficult to implement on a GPU, the calculations were 
performed entirely on a CPU. 
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Figure 3: (a) The difference in free-energy, ^(fc) = &Fknotting{.K) — AFfc MoM / n g(0), against K for 
different tensions, t: 0.\k B T/o (x, black); OAJcbT /o (□, red); 0M B T/o (0> green). Note the 
minimum at fc = K m i n , which decreases for increasing x. (b) The free-energy difference with a term 
proportional to the high A limit in Eq. (2) subtracted: V P( k) — 1 . 1 1 \f%Tt 2 KO% plotted against K for 
the same T. Error bars were estimated by performing three independent repeats of the simulations. 



We first present, in Figure 3(a), results for *¥(k) = AFknotting(K) — AFhwtting(0) a s a function 
of fc for T = 0. 1, 0.4 and 0.8fcgr / a. As expected, we observe that there is a minimum at non-zero 



8 



K, which we denote fc m ,„, and which decreases with increasing tension. In Figure 3(b) we also 
plot the same data subtracting a term proportional to V&n 2 KOX, the expression for AF^notting at 
high A (Eq. (2) with A = ok). The additional proportionality factor of 1.11 was determined by 
fitting AFknotting (k) — AF^ offi „ g (0) for x = 0.4 and O.Sk B T/o for K > \5ksT . For both, the same 
factor was found to the accuracy that is given. The extra factor is likely necessary because our 
polymers do not have negligible thickness. To within errors, the curves for X = 0.4 and O.SksT/a, 
with the expression subtracted, become flat for higher K. This suggests that for these K values 
we have reached the regime where AFknotting is dominated by the bending and tension terms. We 
further observe that, at the position of the minimum of the knotting free energy cost, the quantity 
AFknotting ( K ) ~ AFknottingity — 1.11 V&7t 2 KOX still has a relatively steep slope, confirming that the 
entropic contribution is important in determining the position of the minimum. 
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Figure 4: The minimum value K m j n of bFknatting against the applied tension x for N = 256 (x, 
black) and N = 512 (Q>> red). The solid line is a fit to the five data points for N = 256 with highest 
x values, it has a slope of —0.50 ±0.01. Errorbars were estimated by performing three independent 
repeats of the simulations. 

In Figure 4, we show the dependence of K min on x for N = 256. Plotting on a logarithmic scale, 
we see that the points for the highest five X show a power-law relationship. Fitting to these data, we 
find an exponent of— 0.50±0.01. We thus obtain a power-law dependence of the optimal rigidity 
on the tension that we anticipated in Eq. (1), but with an exponent different than the t = —0.43 we 
found from Figure 1, as expected. For the lowest two x we see that the curve deviates from this 
power-law relationship. This may be attributed to finite size effects. To verify this we repeated 
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simulations for the three lowest % for N = 512: the results are also plotted in Figure 4. We observe 
that, as expected, the results are consistent with the same power-law relationship and also follow it 
to lower x. 

We expect K m j n to be approximately that value of bending rigidity for which the size of the knot 
is equal to the bending length. We consider the variation of the number of the beads in the knot at 
Kmin, Nknot ( K min) , with T. We take K m i n to be given by the best fit relationship from Figure 4. We 
plot the results for A^„ of (fC m; „) in Figure 5. By fitting, we find an exponent of —0.56 ±0.02, close 
to —0.50 ± 0.01 : indeed, Nk n ot( K min) ~ Knin because the polymer is stiff at the scale of the knot. 
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Figure 5: The number of bead in the knot at K m ,„, Nk n ot( K min) against x. The solid line is a fit with 
a slope of —0.56 ± 0.02. Error bars were estimated by performing two independent repeats of the 
simulations. 

We have found that AFknatting has a minimum at a non-zero value of the bending stiffness, 
namely K m [ n . We therefore expect that, if we consider a knotted polymer with non-uniform flexi- 
bility under tension, x, the knot will be more likely to be found in a region with K m i n than in other 
regions. To test this, we consider a polymer of N = 512 beads at x = O.SkgT /o, split into two 
halves: the first 256 beads have Kj = Kq ^ K m i n . The second 256 beads have K",- = 1.806^7' ~ K m i n 
for this x. In Figure 6, we plot results for K"o = 0, 0.4353fcgT ', 0.8706^7 and 3.842fcg7\ i.e. three 
regions with K"o < K m i n and one with Ko > K m i n . Results are binned into 8 bins of 64 beads each. 
In each case we find that the probability of finding the knot in the region with K m i„ is higher. In 
other words, the knot prefers to localize in the region where K w Kq. Furthermore, we find that 
the probabilities are approximately those that would be expected from the free energy calculations. 
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For Kb = in Figure 6, the ratio between the average of the first four bins and that of the sec- 



ond four is 4.9 ±0.5, giving an expected free energy difference of l.6±0.lk B T . The minimum 
AF k nottin g (K) - AF knottlllg (0) for % = 0M B T/a in Figure 3(a) is -1.52±0.02fc 5 7\ 



Figure 6: Probability density, p, of finding the knot at a given position along the polymer under 
tension, % = O.Sk B T/a. For beads 256 — 51 1, K,- = 1.806^7 m K min , whilst for beads — 255 
Ki = (x, black), Kj = 0A353k B T (□ ,red), ?q = O.S106k B T (O, green) or fq = 3M2k B T (A, 
blue). Errorbars were estimated by performing three independent repeats of the simulations. 

To summarize, inspired by correlations between polymer flexibility and knotting seen in biol- 
ogy, we have investigated how the cost of forming a knot in a polymer under tension, t, depends on 
the polymer's stiffness, controlled in our model by K. For high K, our results agree with a simple 
expression including only bending and tension, whilst for lower K entropy must also be taken into 
account. There is a non-zero minimum of the free energy difference between unknotted and knot- 
ted states at k = K m i n . The position of the minimum is seen to depend on tension as fc m ,„ ~ T~ ' 5 . 
We argue that K m i„ is determined by the relative sizes of the knot and the bending length and find 
that the number of polymer beads in the knot at K m j n is consistent with this argument. We consid- 
ered knotted polymers with two sections with different K and found that the knot is more likely to 
be found in the section with K m j n . 

Biological DNA is typically highly confined and in future work it would be interesting to inves- 
tigate the effect of confinement on the results we have observed. 28 ' 29 It would also be interesting to 
investigate how the position of cleavage sites relative to regions of different flexibility affects the 
steady state level of knotting, 30 as well as looking into how the effect of flexibility may combine 
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with previously suggested topoisomerase II guidance mechanisms. 12 Finally, it would be intrigu- 
ing to investigate how non-uniform flexibility affects the diffusional dynamics of a knot along a 
polymer. 31 ' 32 
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